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ABSTRACT 

The direct imaging from the ground of extrasolar planets has become today a major astronomical and biological 
focus. This kind of imaging requires simultaneously the use of a dedicated high performance Adaptive Optics 
[AO] system and a differential imaging camera in order to cancel out the flux coming from the star. In addition, 
the use of sophisticated post-processing techniques is mandatory to achieve the ultimate detection performance 
required. In the framework of the SPHERE project, we present here the development of a new technique, based 
on Maximum A Posteriori [MAP] approach, able to estimate parameters of a faint companion in the vicinity 
of a bright star, using the multi-wavelength images, the AO closed-loop data as well as some knowledge on 
non-common path and differential aberrations. Simulation results show a 10 -5 detectivity at 5a for angular 
separation around 15-^ with only two images. 

Keywords: Image processing, exoplanet detection, differential imaging, inverse problem, regularisation 

1. INTRODUCTION 

Today more than 150 exoplanets have been detected. But a great number among them are known by indirect 
gravitational effects on their parent star. This indirect detection and study allows one to estimate physical 
parameters of the companion, like its orbital period or mass, but does not indicate its atmosphere composition 
or its temperature. Exoplanet direct detection from the ground represents today a great scientific gain on our 
knowledge of exoplanet, since it allows one to perform spectroscopy of the planet. But such a detection needs 
a major improvement of technologies in use, since the star and its companion are separated by a fraction of 
arcsecond, and the flux ratio between them is extremely high (10 6 ). The SPHERE instrument, 1 a VLT Planet 
Finder, will allow to detect photons coming from hot Jupiter planets and will be installed on VLT in 2010. This 
instrument is composed of a high performance extreme AO system, 2 an optimised coronagraphic device 3 and 
a dual band imager. 4 But a dedicated post-processing method is mandatory in order to achieve the ultimate 
detection level of SPHERE. In this paper, we will consider the case of extreme AO coupled to differential imaging. 
The common use of spectral differential images is to perform differences between images at different wavelength in 
order to calibrate the residuals of aberration not corrected by AO and the residuals of diffraction not canceled by 
the coronagraph. The main limitation of differential imaging comes from differential aberrations between the two 
images, or between object images and reference images obtained at different times. The principle of differential 
imaging is detailed in the next section. We propose in the third section an optimised method dedicated to our 
specific issue, based on maximum a posteriori approach and able to estimate the turbulence parameters and the 
object in a pair of images. We present in the fourth section simulation results for the estimation of the turbulent 
phase structure function and the object. 

2. SPECTRAL DIFFERENTIAL IMAGING 

Spectral differential imaging is an instrumental method that aims at "attenuating" the flux of the central star 
with respect to the flux of the potential companion. This method was first initiated by Racine 5 and Marois. 4 
Thus, differential imaging plays a role slightly similar to a coronagraph. The difference between differential 
imaging and a coronagraph is that a coronagraph subtracts only the coherent light to the signal, but before 
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detection. Therefore photon noise in coronagraphic images is also attenuated. The differential imaging is able 
to subtract also the star light, but after detection. The photon noise variance in differential imaging is therefore 
doubled in the combined images. Spectral differential imaging consists in acquiring two simultaneous images of 
a system star-companion at different wavelengths. These two images are rescaled spatially and in intensity and 
combined in a subtraction that reduces the flux of central star and of the residual speckle. Here we will only 
treat differential imaging without coronagraph, for simplicity. 

The subtraction should reduce the star light, but not the companion contribution in the image. This is 
possible if there are strong features in the companion spectrum. In the case of the giant gaseous exoplanets 
searched by the SPHERE project, a strong absorption band due to methane exists at 1.62/xm and can be used 
in such a subtraction: the imaging wavelengths have to be chosen inside and outside the methane band (for 
example 1.575 and 1.625 fJ,m), so that the companion emits at one wavelength, but is drastically less visible at 
the other one. On the other hand, the wavelengths have to be close enough to ensure the speckle pattern of the 
central star only differ in the two images by a spatial and intensity scaling. 

Let us study the image formation theory and the limitation of differential imaging. The expression of the two 
images i\ 1 and i\ 2 can be written as a convolution of the observed object by the Point Spread Function [PSF] 
of the instrument plus additive noises due to photon statistics and electronics: 



i Xl (a) = h Xl (a) * o\ 1 (a) + m 

ix 2 (a) = hx 2 (a) * ox 2 (a) + n 2 (1) 

with h\ 1 and h\ 2 the PSF's in the two imaging channels which depends on turbulence parameters and static 
aberrations in the imaging path, a the angular position in the image field or the object field, o the observed object, 
* stands for the convolution process, n\ and n 2 stand for the noise in the images. For unresolved planets and 
star, the observed object is the sum of Dirac functions weighted by the total flux of the star and the companions 
at their respective position. The images are centred on the star, whose flux is supposed to be the same at the 
two wavelengths. 

o Al (a) = F Q S(a ) + F hi S(di) 

i 

o X2 (a) = F Q S(d ) + ]T F 2ti S(Si) (2) 

i 

with F and a the total flux and position of the central star, F\^, F 2 ^ and on the flux of the different 
companions at wavelengths Ai and A2 and their positions. 

There are two different ways to perform image subtraction: the Single Difference [SD] cancels the effect of 
the common static aberrations, the Double Difference [DD] cancels the effect of both common and differential 
aberrations and is therefore photon noise limited. 

- The SD consist in directly subtracting the two images and allows therefore to cancel the effect of the common 
aberrations in the two imaging channels. The images have to be spatially rescaled at the same wavelength by a 
^ dilation in the focal plane of the second image. 

isD(a) = - U 2 (t^<5) (3) 

M 

If we consider the case where there is no differential aberration, and if Ai and A2 are sufficiently close then 
the PSF h\ 1 can be well approximated by h\ 2 rescaled at Ai. The limitation of this rescaling is Therefore 
the SD gives a good approximation of the difference of the companions convolved by first PSF h\ 1 , as the star 
light has been totally reduced: 



isD(a) = ^F lti 6(ai) - ^ F 2li S(^di)^ * h Xl 



+ ni - n 2 (4) 



But in a more realistic case, the static differential aberrations are not null and the difference between the 
image i\ 1 and the image i\ 2 rescaled at Ai makes appear the effect of differential aberrations. 

The two images have to be acquired simultaneously in order to see the same acquisition conditions (turbulence 
parameters, guide star magnitude, AO performance...). Two imaging channels are therefore used, each of them 
acquiring an image centred on the imaging wavelength. The efficiency of this subtraction depends on differential 
aberration amplitude between the two optical imaging channels, since these aberrations are the main difference 
between the two combined images. 6 

- The DD aims at solving the SD limitation by using two reference images obtained on a reference star with 
the same imaging tool but at another time, the DD therefore cancels the effect of the differential aberrations 
(assuming that they have not evolved between the two observations): 

iDD(a) = (ix^a) -ix 2 - (w,Ai (") - W,a 2 (^f")) ^ 

The reference images are acquired at a different time, and on a different position on sky. This method is 
therefore sensitive to the evolution of the observing conditions between the acquisition of the two pairs of scientific 
and reference images. The evolution of turbulence parameters, AO performance, and most of all the evolution 
of quasi-static aberrations are the main limitations of the DD method. 



3. POST PROCESSING FOR DIFFERENTIAL IMAGING 

As explained before, the detection of low flux companions (contrast around 10 6 between central star and com- 
panion) requires the perfect calibration of both differential static aberrations and system parameters (AO per- 
formance). In a first approximation, we assume that the static aberrations in each imaging channel are perfectly 
known. This is well achieved by using a phase diversity calibration, as described by Sauvage et al.. 7 In this frame- 
work, we present here a new post-processing deconvolution method based on a MAP approach that estimates 
the turbulence-induced PSF and the observed object. 

3.1. Separation static / turbulent aberrations in long exposure images 

The image formation from the ground of stellar objects is perturbed by two factors: the atmospheric turbulence 
and the static aberrations of the telescope. The aberrant pupil phase is therefore the sum of two terms: <j> — (t>t+4>s 
with <j> t the turbulent part and <j) s the static part of the phase. The turbulent phase <j> t is a random variable of 
time and position in pupil plane and is therefore characterised by its structure function whereas the static 
phase <t> 8 does not depend on time and is deterministically known. If the turbulent phase is stationary 8 (as for 
uncorrected turbulence) then it has been shown by Roddier 9 that the OTF is the product of the long exposure 
turbulence-induced OTF and of the static OTF: 

h(f) = cxp (-^D^Xf)^ -L jfjf P(f+ A/) exp (i.4> s {r + A/)) .P(f)*. exp (-i.0 a (A/)) d 2 f (6) 

with 

• P(r) the pupil function 

• D^(Xf) the atmospheric phase structure function after AO correction at wavelength A: 

DtM^dMf+fi-MW (7) 



The phase structure function (A/) is a statistical term that quantifies the turbulent phase variations for 
two points separated by p = A/ in the pupil plane and its shape depends on turbulence parameters and on AO 
performance. If the turbulence is corrected, depends both on f and p 

D^ t (r,p) = (\Mr + p)-Mm (8) 

The average (•) in the expression of D^, t (f,p) is theoretically an average on phase occurrences (and thus on 
time), but may be approximated by an average (■}? on r (stationarity approximation). This simplified expression 
D<p t (p) = (\4>t{r + p) — 4>t{r)\ 2 )r is therefore independent of r. This stationarity approximation is justifiable in 
the case of a Kolmogorov turbulence statistic, and often used also in the case of AO-corrected turbulent phases. 8 

Equation (6) shows that the global OTF is the product of a turbulence-induced OTF and a static OTF: 

Hf) = W)-hs(f) (9) 

with h t (f) the long exposure OTF due to turbulence only, and h s (f) the OTF due to telescope and the 
aberrations. 

The structure function at a wavelength A2 can be rescaled at another wavelength Ai by the operation described 
in Equation 10, in order to compute the turbulence-induced OTF in the first image. Thus, the turbulent OTF 
/it.A x in the first image can be computed thanks to this structure function at A2 by using relation 11. The image 
formation described in Equation 1 can now be rewritten as in Equation 12 taking explicitly into account the 
turbulent and static components of the phase. 

D^xM = (^?D^ tM (^p) (10) 
^A^cxp^-i^) 2 ^,,^^ (11) 



« Al (<S) = /it,Ai («) * /i s ,Ai (a) * o\ 1 (a) 

i\ 2 (a) = h t .x 2 (a) * h s ,\ 2 (a) * o\ 2 (a) (12) 

with ft, ti Ai and h ty \ 2 the turbulent long exposure PSF depending on turbulent phase structure function D ( p i \(p} 1 
h s ,\ 1 and h s ,\ 2 the PSF depending on static aberrations (j) s ,i and <f) s ^ in the two imaging channels. 

3.2. The post-processing framework 

The main limitation of differential imaging comes from differential aberrations in the two spectral channels which 
creates different static pattern in the images in the case of SD, and the evolution of these patterns due to system 
state modification in the case of DD. Therefore in our new approach we propose to estimate the PSF and the 
observed object o in the two images i\ 1 and i\ 2 . The estimation of the PSF reduces to that of its residual 
turbulent component as the static aberrations are supposed to be measured separately. The estimation is done 
thanks to the minimisation of an adequate MAP criterium J(D^,o). 

The MAP approach is based on writing the probability V — P(o, D^i^.i^) of a given object and 
knowing the images using Bayes' theorem (see Equation 13). Finding the best object and structure function 
means maximising the probability V with respect to o and D^. 

V(o,D^) - P(o,D^\i Xl ,i X2 ) cx P(i Xl ,i X2 \o,D^) ■ P{o) ■ P(D^) (13) 

The first factor P(i\ 1 , i\ 2 \o, D^) is called "likelihood term" and embodies the relationship between data and 
the sought parameters. Its statistics is given by the noise statistics in the image (stationary white Gaussian 



noise in a first approximation). The other probabilities are the a priori knowledge we have on the parameters 
to estimate. These regularisation terms allow to smooth the criterium and to accelerate its minimisation. For 
instance, the turbulent phase structure function has a particular shape depending on turbulence parameters and 
may therefore be taken into account in this regularisation term. 

The criterium to minimise is J = — ln('P) with V written in Fourier space and may be rewritten as (Equation 
14). 



+ I \h 2 (/) - h tM (Aa, /) • h sM (/) • 5 A2 (/) || 2 

+ Jr,d^{D 4> ) + Jr, (o) (14) 

where ~ denotes the Fourier transform, Jr^d^ and Jr,o denote regularisation terms accounting for a priori 
knowledge we may have on the parameters to estimate. 

3.3. Assumption and subsequent simplified method 

In the framework of this deconvolution process, we make the following assumption in order to simplify the 
minimisation and demonstrate the feasibility of such a global technique: 

We assume that the companion presents particular spectral signature : it emits light at the first wavelength 
and is totally undetectable at the second wavelength, it means the object is an ideal hot Jupiter and presents 
strong absorption line around 1.6^m. The second image i\ 2 can therefore be seen as a calibration PSF, and the 
global minimisation may be approximated by the three following steps : 

1) Estimation of the structure function Z?0(A 2 • f) in the image i\ 2 without the object, knowing the static 
aberration. This corresponds to the minimisation of the two middle term with respect to in Equation 14: 
likelihood term on i\ 2 , and regularisation term on D^. 

2) Rescaling of the structure function (estimated at A2) at wavelength Ai according to Equation 10 and 
computation of global PSF hi of the first image i\ lt knowing the static aberration of the first channel 
according to Equation 11. 

3) Deconvolution of the first image with the previously computed PSF hi , and estimation of the object in i\ 1 . 
This corresponds to the minimisation of the two terms depending on o only (first and last) in the criterium 
J with respect to o . 

3.4. Estimation of phase structure function 

The turbulent phase structure function gives a statistical knowledge on a turbulent phase. For a turbulent phase 
following a Kolmogorov profile, the structure function is given by the relation D^(p) = (■^■) 3 , with ro the Fried 
parameter. But for a turbulent phase corrected by an AO system, this relation takes into account the AO system 
parameters and is here numerically estimated. The Figure 1 shows typical profiles of D<p for the turbulence and 
AO conditions explained in section 4, with variations of seeing. 

Let us study the first step of the method : the estimation of structure function in a calibration long 
exposure PSF. A criterium (see Equation 15) is used for this minimisation, based on the likelihood term and a 
regularisation term on D<j>. 



JiPt) - INa 2 -F.cxp(--^)-^ A2 || 2 
+ Jd^D*) 



(15) 
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Figure 1. Profiles of phase structure function for a turbulent phase corrected by AO. Condition of simulation : Paranal 
+ SAXO, with different values of seeing. 



with F the flux of the observed star, the structure function to estimate and h Si \ 2 the static PSF due to 
static aberration. is the only estimated parameter, since the star flux as well as the static aberrations (and 
therefore the static PSF) are assumed to be known. 

The regularisation term Jd {D$) is an adaptive smoothness term on estimated designed to avoid noise 
amplification during estimation and to allow the extrapolation of D<p to regions where the static OTF is very 
small (or even null). The regularisation is done using the gradient of structure function VZ?0 and penalises 
deviations between two adjacent pixels according to a typical adaptive variance, depending on pixel position. 
This term is computed as explained in Equation 16. 

Jd, (Aa) - \ (V£>*)* C v \ (VI>,) (16) 

with Cvd^, the covariance matrix of the gradient of phase structure function VD^. This covariance matrix 
has been estimated on different occurences of VD^, these occurences have been generated with different value 
for ro, wind speed or star magnitude. Cyr^ quantifies the typical variability of and allows one to correctly 
weigth the regularisation term. A common issue in regularised inversion methods and criterium minimisation is 
how to choose the hyperparameter, that balances the two terms of the criterium. With the Bayesian apporach 
adapted here and with the use of a Cxjd^ estimated by simulations, there is no such hyper-parameter to be tuned 
and the estimation of is completely unsupervised. 

3.5. Object estimation 

In our procedure, the structure function and the object estimation are done sequentially. This simplified approach 
gives a good idea of the global approach performance, even though a global minimisation should be even more 
precise and therefore lead to a slightly better object estimation (which is the final goal). 



The object estimation is done using MISTRAL 10 algorithm developed at ONERA. This algorithm is based 
on the minimisation of the following criterium, and gives the best object given an image, its PSF and a priori 
knowledge: 

J(o) = ||i Al - h tM * h sM * o|| 2 + J R {o) (17) 

with h t .\ 1 the turbulent PSF at Ai computed with the estimated D^, Jr{o) a regularisation term accounting 
for a priori knowledge on the object. This regularisation term may contains different terms. In our particular 
problematic, we used a positivity constraint and a quadratic linear-quadratic regularisation. 10 

4. RESULTS 

In this section, we validate our post-processing method on simulated data. The simulation conditions are detailed 
in the following list, and correspond to a 8m class Telescope with an Adaptive Optics system of high performance 
like SPHERE, and a turbulence profile corresponding to a typical Paranal sky. The goal of this simulation is to 
compare the detectivity of Single Difference, Double Difference and our approach. 

Conditions : 

• Ai = 1.60/xm, A2 = 1.58/im (corresponding to two wavelengths inside and outside the methane absorption 
line) 

• Turbulence parameter : a typical Cn 2 profile for Paranal is being used with an average wind speed of 12.5 
m/s, and seeing of 0.8 arcsecondes. 

• Adaptive Optics parameters : as extreme- AO. 41x41 actuators with spatially filtered Shack Hartmann 
WFS, a L3CCD working at 1.2kHz sampling frequency. The guide star has a V-magnitude of 8. 

• The static aberration component is randomly generated according to a \ spectrum (n being the radial 
Zcrnike order) with 1300 Zernike coefficients and a differential wavefront error of lOnm RMS in each 
channel, i.e., the total differential wavefront error is 14nm RMS. 

• Imaging parameters : 256x256 images, with a 8m telescope. The different PSF's h\ i are generated at 
Shannon (i.e., one pixel is jfc arcsecond on sky) and are therefore already spatially rescaled. The images 
i\ ± and i\ 2 are generated by convolution of the object and the PSF of each channel. 

• The object at the first wavelength is a star and three companions with a flux ratio of 10~ 3 for the first 

image, and the star alone for the object at the second wavelength. The star flux is set to a total of 10 7 

photons for each wavelength. The companions are located close to the star, respectively at 2.5, 5.0 and 7.5 
x_ 

D ■ 

Such conditions allow us to generate quite realistic images (see example on Figure 2) that are processed by 
our method. 

4.1. estimation : simulation results 

The image i\ 2 is processed in order to estimate the phase structure function via the minimisation of the criterium 
presented in previous section. Figure 3 shows results of D$ estimation. The true (used to generate the images) 
on the left shows the plateau value and central features characteristic of AO system. In the middle, the estimated 
structure function without regularisation (only the likelihood term is used in the criterium). Noise on the edge 
of the circular support of D<p is amplified. The use of adaptive regularisation (on the right) allows us to reduce 
this noise amplification and gives a far better estimation of D^. The error profiles are plotted on Figure 4. 

Without regularisation, the error is lower than 0.03 rad 2 . Figure 4 shows the gain brought by regularisation 
on Dff, estimation : the maximum error for high frequencies is one order of magnitude fainter when regularisation 
is used during estimation. 

The adaptive aspect of the regularisation allows a powerful smoothing of the estimated at the edges, and 
simultaneously a data-driven precise estimation of the quite oscillating D$ near the center. 



Figure 2. Example of spectral images, logarithmic scale, two of the three companions around the central star are visible 
on the left image (Ai = 1.6/im), and only the star in the right image (A2 = 1.58/im). 



4.2. Computation of hi 

The estimated at A2 is now used to compute hi, the PSF of the first image. This computed PSF will then be 
used in the deconvolution of the first image i\ r . Once again we assume to know perfectly the static aberrations 
of first imaging channel. The OTF hi is thus computed as product of turbulent OTF and static OTF by mean 
of Equations (9), (10) and (11). 

4.3. Object estimation 

One can estimate different objects with the different D^, estimated previously : the rough without regulari- 
sation, regularised or by using the true D^, used for images simulation. The different objects estimated are 
gathered in Figure 5, and compared with results of differential imaging in different cases. 

[Top left] The observed object. Central star has a total flux of 10 7 photon, the three companions have a ratio 
of 10 -3 compared to central star and are situated at 2.5, 5.0 and 7.5 inside the AO halo. 

[Top middle] result of single differential imaging. The two images at Ai and A2 are spatially rescaled before 
subtraction. The effect of differential aberrations on central star reduces contrast around it, the first companion 
is unseeable and the second one is visible. 

[Top right] result of double differential imaging. A difference of reference images has been subtracted to the 
single difference of images. This combination of images plays the role of a calibration of differential aberrations 
and allows to reduce their effect. This result is ideal since it does not account for slow variations of static 
aberrations between the two images, of evolution of turbulence parameters between the acquisition of object 
images and reference images. It is therefore a perfect DD, only limited by photon noise. 

[Bottom left] object estimated after deconvolution by the PSF /ii.zv This PSF has been computed with 
estimated D$ with regularisation. The two farest companions are clearly visible, the flux ratio is almost respected. 
The closest companion is quite visible, but a bit hidden by residual flux coming from the central star. 



Figure 3. True D^, estimated without regularisation, estimated with regularisation. 
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Figure 4. X cross-section for [left] true and estimated structure function with and without regularisation, and [right] 
absolute error on estimated structure function with and withoutu regularisation. 



[Bottom right] object estimated after deconvolution with the PSF obtained with regularised D$. The noise 
in this object is slightly fainter than in the previous one, and the central star and the first object are better 
defined. 

5. PERFORMANCE EVALUATION 

The method performance is evaluated here by the 5a detectivity profiles. Detectivity profile is obtained as 5 
times the standard deviation computed azimutally on the result of Single Difference, Double Difference, and 
object estimation by our approach normalised to peak flux in image. 



Figure 5. logarithmic scale of [Top left] observed object, and result of [Top middle] SD and [Top right] DD, [Bottom left] 
deconvolved object with the PSF computed with non-regularised structure functions, [Bottom right] deconvolved object 
with the PSF computed with regularised structure function. 



These detectivity profiles are shown on Figure 6. The simulation conditions are the one listed on section 4. 
As the static differential aberrations are weak (lOnm RMS on each imaging channel), the DD is better than SD 
only close to optical axis (closer than 5-^) when differential aberrations effects dominates. Far from the optical 
axis, the two differences arc photon noise limited and the SD gives therefore slightly better detectivity. We found 
the expected gain in \pi. Whatever the angular separation, our method gives better detectivity. The gain is 
more than 5 in the whole field of view. It is due both to the concentration of the object light in one pixel and 
to photon noise reduction due to our regularised deconvolution. 
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Figure 6. Averaged detection profiles at 5cr in the case of rough image, Single difference, double difference and our 
simplified object estimation. 



6. CONCLUSION 

We propose a method which allows to solve SD limitations (differential aberrations) and gives better results than 
DD, without reference images. In a perfect case, detectivity at 5a reaches less than 10 -5 at 15-^. This result has 
still to be tested in more complex cases (slowly evolving static aberrations or mis-calibration, residual background) 
but gives a rough idea of the potentiallity of the method. The perspectives for this method are to perform a global 
estimation of the parameters (structure function, object in the two images and static aberrations), to process 
real images obtained on a differential imager like NACO SDI, and to generalise our approach to coronagraphy. 
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